BLOCK DATA FRONTF 
include "common.h" 
END 

PROGRAM XSTREAM 
include "common.h" 
REAL*8 X(31,31,15) ,Y(31,31,15) ,Z(31,31,15),GPJ(31,31,15) 
REAL*8 GQJ(31,31,15) ,GRJ(31,31,15) 
CALL UGRID(X, Y, Z) 
CALL GRDGEN(X,Y,Z,GQJ,GPJ,GRJ) 
END 
C end of XSTREAM 

SUBROUTINE CONCAT(RESULT,STRl,STR2) 
CHARACTER *(*) RESULT,STR1,STR2 
INTEGER IPOS,IP2 
C 
C SEARCH STRl FOR FIRST EMPTY SPACE 
C 
IPOS = 1
IP2 = 1
100 CONTINUE 
RESULT(IP2:IP2)=STR1(IPOS:IPOS)
IPOS = IPOS + 1
IP2 = IP2 + 1
IF((STRl(IPOS:IPOS+l).NE.'').AND.(IPOS.LE.LEN(STRl))) GOT0 100 

C 
C CONCATENATE THE SECOND STRING
C
IPOS = 1 
200 CONTINUE
RESULT(IP2:IP2)=STR2(IPOS:IPOS)
IPOS = IPOS + 1
IP2 = IP2 +1
IF(IPOS.LE.LEN(STR2)) GOT0 200 

C 
C PAD WITH SPACES 
C
300 CONTINUE
RESULT(IP2:IP2)=' '
IP2=IP2 + 1
IF(IP2.LE.LEN(RESULT)) GOTO 300
C 
RETURN 
END 


SUBROUTINE SAVEGD(X,Y,Z)
include "common.h"

C       This subroutine writes the geometric data into a file for later use
C       The file called fname.grd

REAL*8 X(NIP1,NJP1,NKP1),Y(NIP1,NJP1,NKP1),Z(NIP1,NJP1,NKP1)
CHARACTER*50 FNAME

C OPEN THE FILE C 
CALL CONCAT(FNAME,FILNAM, '.out') 
OPEN(22,FILE=FNAME,FORM='UNFORMATTED') 

C NOW WRITE THE FILE HEADER 
C 
WRITE(22) NI,NJ,NK 
C 
C AND THE GRID POINTS 
C 
WRITE(22) X 
WRITE(22) Y 
WRITE(22) Z 
C
CLOSE(22)
C 
RETURN 
END 


SUBROUTINE GRDGEN(X,Y,Z,GPJ,GQJ,GRJ) 
include "common.h" 
INTEGER WHAT, LOOP, I, J, K 
REAL*8 DS,RTDS,TDS,X(NIP1,NIP1,NXP1),Y(NIP1,NIP1,NKP1)
REAL*8 Z(NIP1,NIP1,NKP1),GQJ(NIP1,NIP1,NKP1),GPJ(NIP1,NIP1,NKP1)
REAL*8 GRJ(NIP1,NJP1,NKP1)
REAL*8 accn

LOOP = 40 
accn = 0.4D0 
DS=0.1D0 
RTDS = 0.0D0 

C       MAIN LOOP STARTS HERE ... 
1       CONTINUE 

C 	CALL SAVEGD(X,Y,Z) CALL SAVEGRD (X ,Y,Z) CALL SAVEMTV(X,Y,Z) 
PRINT*,'**** GRID GENERATOR MENU ****' 
PRINT*, '1 Change number of iterations (currently ',LOOP,' ) ' 
PRINT*,'2 Change ds (wall distance -currently ',DS, ') ' 
PRINT*,'3 Gradually decreasing ds '
PRINT*,'4 Elliptic generator + boundary clustering' 
PRINT*,'5 Change acceleration parameter' 
PRINT*,'6 Quit' 
READ*,WHAT 

GOT0 (10,20,30,40,50,60) WHAT 
GOT0 1 

10 	CONTINUE
PRINT*, 'Enter no. of iterations ... '
READ*,LOOP 
GOT0 1

20 	CONTINUE
PRINT*, 'Enter ds ...' 
READ*, DS
TDS = 0.0D0
RTDS = 0.0D0
GOT0 1 

30 	CONTINUE
PRINT*, 'Enter target DS ..'
READ*,TDS
PRINT*, 'Enter fractional decrease in ds aer step ..' 

40 	CONTINUE 
CALL ELLIPT (X, Y, Z, LOOP, 1.0D0, GPJ, GQJ, GRJ, DS, accn, RTDS,TDS) 
GOT0 1 

50 	CONTINUE 
PRINT*, 'Enter accn. '
READ*, accn 
GOT0 1 

60 	CONTINUE
RETURN 
END 

SUBROUTINE ELLIPT (X,Y,Z, LOOP, BC, GPJ, GQJ, GRJ, DS, accn, B.TDS, TDS)
include "common. h"
REAL*8 X(NIP1,NJP1,NKP1),Y(NIP1,NJP1,NKP1),Z(NIP1,NJP1,NKP1)
REAL*8 BC,accn,xn,yn,zn,xe,ye,ze,xz,yz,zz,xen,yen,zen 
REAL*8 xez,yez,zez,xnz,ynz,znz
REAL*8 gll,g12,g13,g21,g22,g23,g31,g32,g33
REAL*8 GQJ(NIP1,NJP1,NKP1),GPJ(NIP1,NJP1,NKP1),GRJ(NIP1,NJP1,NKP1) 
REAL*8 dummyl, dummy2, dmy3, gapl, Gp, Gq, Gr, IDS, RTDS, TDS, TZRES 
REAL*8 xAterm, yAterm, zAterm, accnx, accny, accnz, DS, XRES, YRES, ZRES 
REAL*8 JAC,a11,a12,a13,a23,a33,a22,sum,dd,dx,dy,dz,TXRES,TYRES 
REAL*8 AX(100),AY(100),AZ(100),BX(100) BY(100),BZ(100)
REAL*8 CX(100),CY(100),CZ(100),RX(100),RY(100),RZ(100) 
DOUBLE PRECISION cpu1,cpu2,cputime 
INTEGER LOOP,L,I,J,K,ITER 
CHARACTER*l STAT 

STAT=' N' 
accnx = accn 
accny = accn 
accnz = accn 
TXRES = 1.0D0 
TYRES = 1.0D0 
TZRES = 1.0D0 
IDS=DS 

cpul=X05BAF() 

DO L=1,LOOP 
IF ((ABS(TZRES).LT.0.01D0).AND.(RTDS.NE.0.D0)) THEN
IF (DS.GT.TDS) THEN 
DS = DS+RTDS*(TDS-DS) 
ENDIF 
ENDIF 

IF (BC. GT. 0.0D0) THEN 
CALL GRDPQ (X, Y,Z,GPJ,GQJ,GRJ,DS,L) 
CALL SPRDPQ(GPJ,GQJ,GRJ) 
ENDIF 

XRES=0.D0
YRES=0.D0
ZRES=0.D0
TXRES=0.D0
TYRES=0.D0
TZRES=0.D0

C       K SWEEP from lower wall to upper wall
DO I=2,NI
DO J=2,NJ
DO K=3,NK-1 

xn = (X(I,J+1,K)-X(I,J-1,K))*0.5D0
yn = (Y(I,J+1,K)-Y(I,J-1,K))*0.5D0
zn = (Z(I,J+1,K)-Z(I,J-1,K))*0.5D0

xe = (X(I+1,J,K)-X(I-1,J,K))*0.5D0
ye = (Y(I+1,J,K)-Y(I-1,J,K))*0.5D0
ze = (Z(I+1,J,K)-Z(I-1,J,K))*0.5D0

xz = (X(I,J,K+1)-X(I,J,K-1))*0.5D0
yz = (Y(I,J,K+1)-Y(I,J,K-1))*0.5D0
zz = (Z(I,J,K+1)-Z(I,J,K-1))*0.5D0

xen=(X(I+1,J+1,K)-X(I-1,J+1,K)+X(I-1,J-1,K)-X(I+1,J-1,K))*0.25D0
yen=(Y(I+1,J+1,K)-Y(I-1,J+1,K)+Y(I-1,J-1,K)-Y(I+1,J-1,K))*0.25D0
zen=(Z(I+1,J+1,K)-Z(I-1,J+1,K)+Z(I-1,J-1,K)-Z(I+1,J-1,K))*0.25D0

xez=(X(I+1,J,K+1)-X(I-1,J,K+1)+X(I-1,J,K-1)-X(I+1,J,K-1))*0.25D0
yez=(Y(I+1,J,K+1)-Y(I-1,J,K+1)+Y(I-1,J,K-1)-Y(I+1,J,K-1))*0.25D0
zez=(Z(I+1,J,K+1)-Z(I-1,J,K+1)+Z(I-1,J,K-1)-Z(I+1,J,K-1))*0.25D0

xnz=(X(I,J+1,K+1)-X(I,J-1,K+1)+X(I,J-1,K-1)-X(I,J+1,K-1))*0.25D0
ynz=(Y(I,J+1,K+1)-Y(I,J-1,K+1)+Y(I,J-1,K-1)-Y(I,J+1,K-1))*0.25D0
znz=(Z(I,J+1,K+1)-Z(I,J-1,K+1)+Z(I,J-1,K-1)-Z(I,J+1,K-1))*0.25D0

gll = yn*zz-yz*zn
g12 = yz*ze-ye*zz
g13 = ye*zn-yn*ze

g2l = xn*zn-xz*yn
g23 = xn*ze-xe*zn
g22 = xe*zz-xz*ze

g3l = xn*yz-xz*yn
g32 = xz*ye-xe*yz
g33 = xe*yn-xn*ye

a11 = g11*g11 + g21*g21 + g31*g31
a22 = g12*g12 + g22*g22 + g32*g32
a33 = g13*g13 + g23*g23 + g33*g33

a12 = g11*g12 + g21*g22 + g3l*g32
a13 = g11*g13 + g21*g23 + g31tg33
a23 = g12*g13 + g22*g23 + g32*g33

JAC = xe*(gll)-xn*(-g12)+xz*(g13)

Gp=GPJ(I,J,K)
Gq=GQJ(I,J,K)
Gr=GRJ(I,J,K)

xAterm = -(JAC*JAC)*(Gp*xe+Gq*xn+GR*xz)
xAterm = -(JAC*JAC)*(Gp*ye+Gq*yn+GR*yz)
xAterm = -(JAC*JAC)*(Gp*ze+Gq*zn+GR*zz)

C --for x 
dummyl=-a11*(X(I+l,J,K)+X(I-1,J,K)) 
dummy2=-a22*(X(I,J+l,K)+X(I,J-l,K)) 
gapl=-2.0D0*(al2*xen+al3*xez+a23*xnz)


C ---Assemble TDMA Coefficients for x
AX(K)=a33
BX(K)=-2.0D0*(a11+a22+a33)
CX(K)=a33
RX(K)=dummyl+dummy2+gapl+xAterm

C ---for y 
dummyl=-a11*(Y(I+1,J,K)+Y(I-1,J,K))
dummy2=-a22*(Y(I,J+1,K)+Y(I,J-1,K))
gapl=-2.0D0*(al2*yen+a13*yez+a23*ynz)

C ---Assemble TDMA Coefficients for y
AY(K)=AX(K)
BY(K)=BX(K)
CY(K)=CX(K)
RY(K)=dummyl+dummy2+gapl+yAterm

C ---for z 
dummyl=-a11*(Z (I+1, J,K)+Z (1-1, J,K)) 
dummy2=-a22*(Z(I,J+1,K)+Z(I,J-1,K)) 
gapl=-2.0D0*(al2*zen+a13*zez+a23*znz) 

C ---Assemble TDMA Coefficients for z

AZ(K)=AX(K)
BZ(K)=BX(K)
CZ(K)=CX(K)
RZ(K)=dummyl+dummy2+gapl+zAterm 

ENDDO 
C --- end of J loop 

C Implicit BCs for X
C lower wall

RX(3)= RX(3)-AX(3)*X(I,J,2)
AX(3)= 0.0D0

C upper wall
RX(NK-1) = RX(NK-1)-CX(NK-1)*X(I,J,NK)
CX (NK-1) = 0.0D0 

C Implicit BCs for Y 
C       lower wall 

C upper wall

RY(NK-1) = RY(NK-1)-CY(NK-l)*Y(I,J,NK)
CY(NK-1) = 0.0D0

C       --Implicit BCs for Z
C       --lower wall

RZ(3) = RZ(3)-AZ(3) *Z(I,J,2) AZ(3) = 0.0D0 

C --SOLVE FOR X 
CALL TDMA(AX,BX,CX,RX,3,NK-1,100) 

C --SOLVE FOR Y 
CALL TDMA(AY,BY,CY,RY,3,NK-1,100) 

C --SOLVE FOR Z 
CALL TDMA(AZ,BZ,CZ,RZ,3,NK-1,100) 


C       ---Update x,y and z 

DO K=3,NK-1

Z(I,J,K)= accnz*Z(I,J,K)+(1.0D0-accnz)*RZ(K)
if ((j.eq.2).or.(j.eq.NJ).and.(i.eq.2).or.(i.eq.NI)) then 
C -- do nothing 
else if ((j.ne.NJ).and.(j.ne.2)) then
Y(I,J,K)=accny*Y(I,J,K)+(1.0D0-accny)*RY(K)
endif 

if((j.eq.2).or.(j.eq.NJ).and.(i.eq.2).or.(i.eq.NI)) then 
C -- do nothing
else if ((i.ne.2).and.(i.ne.NI)) then
        X(I,J,K)=accnx*X(I,J,K)+(1.0D0-accnx)*RX(K)
endif

ENDDO 
ENDDO 
ENDDO 

PRINT*,' ITER =',L 
CALL GRDEXT(X,Y,Z) 
IF(L.eq.LOOP) THEN 

C --Residual Calculation
DO K=3,NK-1
DO J=3,NJ-1
DO I=3,NI-1 

xn = (X(I,J+1,K)-X(I,J-1,K))*0.5D0
yn = (Y(I,J+1,K)-Y(I,J-1,K))*0.5D0
zn = (Z(I,J+1,K)-Z(I,J-1,K))*0.5D0

xe = (X(I+1,J,K)-X(I-1,J,K))*0.5D0
ye = (Y(I+1,J,K)-Y(I-1,J,K))*0.5D0
ze = (Z(I+1,J,K)-Z(I-1,J,K))*0.5D0

xz = (X(I,J,K+1)-X(I,J,K-1))*0.5D0
yz = (Y(I,J,K+1)-Y(I,J,K-1))*0.5D0
zz = (Z(I,J,K+1)-Z(I,J,K-1))*0.5D0

xen=(X(I+1,J+1,K)-X(I-1,J+1,K)+X(I-1,J-1,K)-X(I+1,J-1,K))*0.25D0
yen=(Y(I+1,J+1,K)-Y(I-1,J+1,K)+Y(I-1,J-1,K)-Y(I+1,J-1,K))*0.25D0
zen=(Z(I+1,J+1,K)-Z(I-1,J+1,K)+Z(I-1,J-1,K)-Z(I+1,J-1,K))*0.25D0

xez=(X(I+1,J,K+1)-X(I-1,J,K+1)+X(I-1,J,K-1)-X(I+1,J,K-1))*0.25D0
yez=(Y(I+1,J,K+1)-Y(I-1,J,K+1)+Y(I-1,J,K-1)-Y(I+1,J,K-1))*0.25D0
zez=(Z(I+1,J,K+1)-Z(I-1,J,K+1)+Z(I-1,J,K-1)-Z(I+1,J,K-1))*0.25D0

xnz=(X(I,J+1,K+1)-X(I,J-1,K+1)+X(I,J-1,K-1)-X(I,J+1,K-1))*0.25D0
ynz=(Y(I,J+1,K+1)-Y(I,J-1,K+1)+Y(I,J-1,K-1)-Y(I,J+1,K-1))*0.25D0
znz=(Z(I,J+1,K+1)-Z(I,J-1,K+1)+Z(I,J-1,K-1)-Z(I,J+1,K-1))*0.25D0

gll = yn*zz-yz*zn
g12 = yz*ze-ye*zz
g13 = ye*zn-yn*ze

g2l = xn*zn-xn*zz
g23 = xn*ze-xe*zn
g22 = xe*zz-xz*ze

g3l = xn*yz-xz*yn
g32 = xz*ye-xe*yz
g33 = xe*yn-xn*ye

a11 = g11*g11 + g21*g21 + g31*g31
a22 = g12*g12 + g22*g22 + g32*g32
a33 = g13*g13 + g23*g23 + g33*g33

a12 = g11*g12 + g21*g22 + g3l*g32
a13 = g11*g13 + g21*g23 + g31tg33
a23 = g12*g13 + g22*g23 + g32*g33

JAC = xe*(gll)-xn*(-g12)+xz*(g13)

Gp=GPJ(I,J,K)
Gq=GQJ(I,J,K)
Gr=GRJ(I,J,K)

xAterm = (JAC*JAC)*(Gp*xe+Gq*xn+GR*xz)
xAterm = (JAC*JAC)*(Gp*ye+Gq*yn+GR*yz)
xAterm = (JAC*JAC)*(Gp*ze+Gq*zn+GR*zz)

C --	residual for x 
dummyl= (X(I+1,J,K)-2.0D0*X(I,J,K)+X(I-1,J,K))*a11 
dummy2=(X(1,J+1,K)-2.0D0*X(I,J,K)+X(I,J-1,K))*a22 
dummy3=(X(I,J,K+1)-2.0D0*X(I,J,K)+X(I,J,K-1))*a33

gapl = 2.0D0*(a12*xen+a13*xez+a23*xnz)
sum=(dummyl+dmy2+dummy3)+(gapl+xAterm)
XRES = XRES + ABS(sum)

C --	residual for y
dummyl=(Y(I+1,J,K)-2.0D0*Y(I,J,K)+Y(I-1,J,K))*a11
dummy2=(Y(I,J+1,K)-2.0D0*Y(I,J,K)+Y(I,J-1,K))*a22
dummy3=(Y(I,J,K+l)-2.0D0*Y(I,J,K)+Y(I,J,K-1))*a33

gapl = 2.0D0*(a12*yen+al3*yez+a23*ynz)
sum=(dummyl+dmy2+dummy3)+(gapl+yAterm)
YRES = YRES + ABS (sum)

C --	residual for z
dummyl=(Z(I+1,J,K)-2.0D0*Z(I,J,K)+Z(I-1,J,K))*a11
dummy2=(Z(I,J+1,K)-2.0D0*Z(I,J,K)+Z(I,J-1,K))*a22
dummy3=(Z(I,J,K+1)-2.0D0*Z(I,J,K)+Z(I,J,K-1))*a33
gapl = 2.0D0*(al2*zen+a13*zez+a23*znz)
sum=(dummyl+dummy2+dummy3)+(gapl+zAterm)
ZRES = ZRES + ABS (sum) 
ENDDO 
ENDDO 
ENDDO 

TXRES =XRES/((NI-1)-2)/((NJ-1)-2)/((NK-1)-2)
TYRES =YRES/((NI-1)-2)/((NJ-1)-2)/((NK-1)-2)
TZRES =ZRES/((NI-1)-2)/((NJ-1)-2)/((NK-1)-2) 

if ((ABS(TXRES).lt.0.01D0).and.(ABS(TYRES).lt.0.01D0) .and. (ABS(TZRES1 .lt.0.01D0) then 
        if (STAT.ne.'C') then
                STAT='C'
                ITER=L 
endif
else
        STAT='N' 
endif  
        PRINT*, TXRES, ' ' ,TYRES,' ',TZRES
endif
PRINT* 


ax=X(24,24,3)-(24,24,2) 
dy=Y(24,24,3)-Y(24,24,2) 
dz=Z(24,24,3)-Z(24,24,2) 
dd=DSQRT(dx*dx+dy*dy+dz*dz) 
ENDDO 

C --	End of Iteration Loop.
cpu2=X05BAF()
cputime=cpu2-cpul
PRINT*, 'CPU time =',cputime
PRINT*,' ITER =' ,ITER 
C --	check spacing ... 

CALL SEESPC (X, Y, Z) 

RETURN 
END

C --	End of Ellipt 




SUBROUTINE TDMA(AA,BB,CC,FR,BL,NL,WL) 
INTEGER L,NL,BL,MAXL 
REAL*8 AA(MAXL),BB(WL),CC(ML),RR(MAXL),TEMP 

C 	Forward elimination 

DO L=BL+1,NL 
TEMP = AA(L) /BB(L-1) 
BB (L)=BB (L) -CC (L-1) *TEMP 
RR(L)=RR(L)-RR(L-1) *TEMP 
ENDDO

C Back substitution
RR(NL) = RR(NL) /BB(NL) 
DO L=NL-1,BL,-1 
RR(L)=(RR(L)-CC(L)*RR(L+1)/BB(L) 
ENDDO 
RETURN 
END 


SUBROUTINE GRDPQ(X,Y,Z,GPJ,GQJ,GRJ,DS,L) 
inc1ude "common.h" 
REAL*8 X(NIP1,NJP1,NKP1),Y(NIP1,NJP1,NKP1),GQJ(NIP1,NJP1,NKP1),DS
REAL*8 GPJ(NIP1,NJP1,NKP1),Z(NIP1,NIP1,NKP1),GRJ(NIP1,NJP1,NKP1)
INTEGER I,J,K,L 

K=2
DO J=2,NJ
DO I=2,NI 
CALL SPQR(GPJ,GQJ,GRJ,X,Y,Z,I,J,K,1,L,DS) 
ENDDO 
ENDDO 
RETURN 
END 

SUBROUTINE SPQR(GPJ,GQJ,GRJ,X,Y,Z,I,J,K,PJ,L,DS) 
include "common.h" 
INTEGER I,J,K,PJ,L 
REAL*8 GPJ(NIP1,NJP1,NKP1),GQJ(NIP1,NJP1,NKP1),X(NIP1,NJP1,NKP1),Y(NIP1,NIP1,NKP1)
REAL*8 Z(NIP1,NJP1,NKP1),GRJ(NIP1,NJP1,NKP1)
REAL*8 xe,ye,ze,xn,yn,zn,xz,yz,zz,DS,sgn,JAC,tmp
REAL*8 xee,yee,zee,xnn,ynn,znn,xzz,yzz,zzz
REAL*8 tl,t2,t3,xen,yen,zen,xez,yez,zez,xnz,ynz,znz
REAL*8 rxe,rye,rze,rxn,ryn,rzn,lxe,lye,lze,lxn,lyn,lzn
REAL*8 rxz,ryz,rzz,lxz,lyz,lzz,xRHS,yRHS,zRHS
REAL*8 gll,g12,gl3,g21,g22,g23,g31,g32,g33
REAL*8 rg13,rg23,rg33,lg13,lg23,lg33
REAL*8 a11,a22,a33,a12,a13,a23 

sgn = DBLE(PJ)

xn = (X(I,J+1,K)-X(I,J-1,K))*0.5D0
yn = (Y(I,J+1,K)-Y(I,J-1,K))*0.5D0
zn = (Z(I,J+1,K)-Z(I,J-1,K))*0.5D0
xe = (X(I+1,J,K)-X(I-1,J,K))*0.5D0
ye = (Y(I+1,J,K)-Y(I-1,J,K))*0.5D0
ze = (Z(I+1,J,K)-Z(I-1,J,K))*0.5D0
g13 = ye*zn-yn*ze
g33 = xe*yn-xn*ye
g23 = xn*ze-xe*zn

tmp = 1.0D0/g33

zz = DS/DSQRT((g13*tmp)*(g13*tmp)+(g23*tmp)*(g23*tmp)+1.0D0)
xz = g13*zz*tmp
yz = g23*zz*tmp

g11 = yn*zz-yz*zn 
g12 = yz*ze-ye*zz 
g21 = xz*zn-xn*zz 
g22 = xe*zz-xz*ze 
g31 = xn*yz-xz*yn 
g32 = xz*ye-xe*yz 

a11 = gll*gll + g21*g21 + g31*g31
a22 = g12*g12 + g22*g22 + g32*g32
a33 = g13*g13 + g23*g23 + g33*g33
a12 = gll*g12 + g21*g22 + g31*g32
a13 = g11*g13 + g21*g23 + g31*g33 
a23 = g12*g13 + g22*g23 + g32*g33

JAC = 1.0D0/(xe*(gll) -xn*(-g12) +xz*(g13))

C - second derivatives at lower wall
tmp = 8.0D0*X(I,J,K+PJ)-7.0D0*X(I,J,K)-X(I,J,K+2*PJ)
xzz = tmp*0.5D0-sgn*3.0D0*xz

tmp = 8.0D0*Y(I,J,K+PJ)-7.0D0*Y(I,J,K)-Y(I,J,K+2*PJ)
yzz = tmp*0.5D0-sgn*3.0D0*yz

tmp = 8.0D0*Z(I,J,K+PJ)-7.0D0*Z(I,J,K)-Z(1,J,K+2*PJ)
zzz = tmp*0.5D0-sgn*3.0D0*zz

xee=X(I+1,J,K)+X(I-1,J,K)-2.0D0*X(I,J,K) 
yee=Y(I+1,J,K)+Y(I-1,J,K)-2.0D0*Y(I,J,K) 
zee=Z(I+1,J,K)+Z(I-1,J,K)-2.0D0*Z(I,J,K)

xnn=x(I,J+1,K)+X(I,J-1,K)-2.0D0*X(I,J,K) 
ynn=Y(I,J+1,K)+Y(I,J-1,K)-2.0D0*Y(I,J,K) 
znn=Z(I,J+1,K)+Z(I,J-1,K)-2.0D0*Z(I,J,K)


C       - xez cross derivatives 

rxe=X(I+1,J,K)-X(I,J,K) 
rye=Y(I+1,J,K)-Y(I,J,K) 
rze=Z(I+1,J,K)-Z(I,J,K) 
rxn=((X(I+1,J+1,K)-X(I+1,J-1,K))+(X(I,J+1,K)-X(I,J-1,K)))*0.25D0 
ryn=((Y[I+1,J+1,K)-Y(I+1,J-1,K))+(Y(I,J+1,K)-Y(I,J-1,K)))*0.25D0
rzn=((Z(I+1,J+1,K)-Z(I+1,J-1,K))+(Z(I,J+1,K)-Z(I,J-1,K)))*0.25D0 

rg13 = rye*rzn-ryn*rze
rg23 = rxn*rze-rxe*rzn
rg33 = rxe*ryn-rxn*rye

rzz=DS/DSQRT((rg13/rg33)+(rg23/rg33)+1.0D0)
rxz=(rg13/rg33)*rzz
ryz=(rg23/rg33)*rzz

lxe=X(I,J,K)-X(I-1,J,K)
lye=Y(I,J,K)-Y(I-1,J,K)
lze=Z(I,J,K)-Z(I-1,J,K)
lxn=((X(I,J+1,K)-X(I,J-1,K))+(X(I-1,J+1,K)-X(I-1,J-1,K)))*0.25D0 
lyn=((Y[I,J+1,K)-Y(I,J-1,K))+(Y(I-1,J+1,K)-Y(I-1,J-1,K)))*0.25D0
lzn=((Z(I,J+1,K)-Z(I,J-1,K))+(Z(I-1,J+1,K)-Z(I-1,J-1,K)))*0.25D0 

lg13 = lye*lzn -lyn*lze
lg23 = lxn*lze -lxe*lzn
lg33 = lxe*lyn -lxn*lye

lzz=DS/DSQRT((lg13/lg33)+(lg23/lg33)+1.0D0)
lxz=(lg13/lg33)*lzz
lyz=(lg23/lg33)*lzz

xez=rxz-lxz
yez=ryz-lyz
zez=rzz-lzz

C       - nz cross derivatives 

rxn=X(I,J+1,K)-X(I,J,K)
ryn=Y(I,J+1,K)-Y(I,J,K)
rzn=Z(I,J+1,K)-Z(I,J,K)
rxe=((X(I+1,J+1,K)-X(I-1,J+1,K))+(X(I+1,J,K)-X(I-1,J,K)))*0.25D0
rye=((Y(I+1,J+1,K)-Y(I-1,J+1,K))+(Y(I+1,J,K)-Y(I-1,J,K)))*0.25D0
rze=((Z(I+1,J+1,K)-Z(I-1,J+1,K))+(Z(I+l,J,K)-Z(I-1,J,K)))*0.25D0

rg13 = rye*rzn-ryn*rze
rg23 = rxn*rze-rxe*rzn
rg33 = rxe*ryn-rxn*rye

tmp = 1.0D0/rg33

rzz = DS/DSQRT((rg13*tmp)*(rg13*tmp)+(rg23*tmp)*(rg23*tmp)+1.0D0)
rxz = (rg13/rg33)*rzz
ryz = (rg23/rg33)*rzz

lxn=X(I,J,K)-X(I,J-1,K)
lyn=Y(I,J,K)-Y(I,J-1,K)
lzn=Z(I,J,K)-Z(I,J-1,K)
lxe=((X(I+1,J,K)-X(I-1,J,K))+(X(I+1,J-1,K)-X(I-1,J-1,K)))*0.25D0 
lye=((Y(I+1,J,K)-Y(I-1,J,K))+(Y(I+1,J-1,K)-Y(I-1,J-1,K)))*0.25D0
lze=((Z(I+1,J,K)-Z(I-1,J,K))+(Z(I+1,J-1,K)-Z(I-1,J-1,K)))*0.25D0 

lg13 = lye*lzn -lyn*lze
lg23 = lxn*lze -lxe*lzn
lg33 = lxe*lyn -lxn*lye

tmp = 1.0D0/rg33

lzz = DS/DSQRT((lg13*tmp)*(lg13*tmp)+(lg23*tmp)*(lg23*tmp)+1.0D0)
rxz = (lg13/lg33)*lzz
ryz = (lg23/lg33)*lzz

xnz=rxz-lxz
ynz=ryz-lyz
znz=rzz-lzz


C en cross derivative
xen=(X(I+1,J+1,K)-X(I-1,J+1,K)+X(I-1,J-1,K)-X(I+l,J-1,K)) *0.25D0
xen=(Y(I+1,J+1,K)-Y(I-1,J+1,K)+Y(I-1,J-1,K)-Y(I+l,J-1,K)) *0.25D0
xen=(Y(I+1,J+1,K)-Z(I-1,J+1,K)+Z(I-1,J-1,K)-Z(I+l,J-1,K)) *0.25D0

C P,Q,R value from RHS times Transpose of the cofactor matrix 
C divided by Jacobian 

xRHS=-JAC*JAC*(a11*xee*a22*xnn+a33*xzz+2.0D0*(a12*xen+a13*xez+a23*xnz)
yRHS=-JAC*JAC*(a11*yee*a22*ynn+a33*yzz+2.0D0*(a12*yen+a13*yez+a23*ynz)
zRHS=-JAC*JAC*(a11*zee*a22*znn+a33*zzz+2.0D0*(a12*zen+a13*zez+a23*znz)


t1 = (xRHS*g11+yRHS*g21+zRHS*g31)*JAC
t2 = (xRHS*g12+yRHS*g22+zRHS*g32)*JAC
t3 = (xRHS*g13+yRHS*g23+zRHS*g33)*JAC

C Updating the P, Q, and R values 

GPJ(I,J,K)=GPJ(I,J,K)+0.05D0*(t1-GPJ(I,J,K))
GQJ(I,J,K)=GQJ(I,J,K)+0.05D0*(t1-GQJ(I,J,K))
GRJ(I,J,K)=GRJ(I,J,K)+0.05D0*(t1-GRJ(I,J,K))



C PRINT*, 'Gp =',GPJ(I,J,K)
C PRINT*, 'Gq =',GQJ(I,J,K) 
C PRINT*, 'Gr =',GRJ(I,J,K) 
C READ* 

RETURN 
END 

SUBROUTINE SPRDPQ(GPJ,GQJ,GRJ) 
include "common.h" 
REAL*8 GPJ(NIP1,NJP1,NKP1),GQJ(NIP1,NJP1,NKP1),GRJ(NIP1,NJP1,NKP1)
INTEGER I,J,K 

DO K=3,NK 
DO J=2,NJ 
DO I=2,NI 
GPJ(I,J,K)=GPJ(I,J,2)*EXP(-0.5D0*(K-2)) 
GQJ(I,J,K)=GQJ(I,J,2)*EXP(-0.5D0*(K-2)) 
GRJ(I,J,K)=GRJ(I,J,2)*EXP(-0.5D0*(K-2)) 
ENDDO 
ENDDO 
ENDDO 
RETURN 
END 

REAL*8 FUNCTION ZCURV(XP,YP) 
REAL*8 XP,YP,CURV 
CURV=l.0D0*EXP(-0.5D0*(XP*XP+YP*YP)) 
ZCURV=CURV 
RETURN 
END 

SUBROUTINE UGRID(X,Y,Z) 
include "common.h" 
C 
C THIS SUBROUTINE IS USED TO SPECIFY THE GRID TO BE USED FOR THE 
C COMPUTATION. EACH OF THE COORDINATES MUST BE SPECIFIED FOR THE 
C GRID BOUNDING THE CONTROL VOLUMES. 

REAL*8 XMIN,XMAX,YMIN,YMAX,ZMIN,ZMAX,XO,XP,XR,PW,RI1,RI2,RJl,RJ2 
REAL*8 X1(NIP1,NJP1,NKP1),Yl(NIP1,NIP1,NKP1) 
REAL*8 X2(NIP1,NIP1,NKP1),Y2(NIP1,NJP1,NKP1),DX1,TSUM,SUM,DELTAY,Yo,Yp,DELTAZ 
REAL*8 Z1(NIP1,NJP1,NKP1),Z2(NIP1,NIP1,NKP1) 

INTEGER CN,CI,CJ 
INTEGER I,J,K 
PARAMETER (XMIN = -4.0D0, XMAX = 4.0D0) 
PARAMETER (YMIN = -4.0D0, YMAX = 4.0D0) 
PARAMETER (ZMIN = 0.0D0, ZMAX=4.0D0) 
DELTAX = (XMAX-XMIN)/DBLE(NI-2) 
DELTAY = (YMAX-YMIN)/DBLE(NJ-2) 
DELTAZ = (ZMAX-ZMIN)/DBLE(NK-2) 

C ** Bottom boundary points redistribution to attract points to the centre of the bump using
C -- geometrical ratio XR ****

CN=INT((NI-1)/2)+1
XR=0.95D0

SUM=0.0D0 
DO I=1,CN-2 
CALL POW(XR,I,PW) 
SUM=SUM+(1.0D0/PW) 
PW=0.0D0 
ENDDO 
        TSUM=2.0D0*(SUM+1.0D0) 
        DX1=(XMAX-XMIN)/TSUM 
        DX(CN)=DX1 
        DX(CN-1)=DX1 

DO I=1,CN-2 
CALL POW(XR,I,PW) 
DX(CN+I)=DX1/PW 
PW=0.0D0 
ENDDO

DO I=CN-2,1,-1 
CALL POW(XR,I,PW) 
DX( (CN-1)-I)=DX1/PW 
PW=0.0D0 
ENDDO 

K=2 
DO I=2,NI 
Y(I,2,K)=YMIN 
Y(I,NJ,K)=YMAX 
ENDDO 

K=2 
DO J=2,NJ 
X(2,J,K)=XMIN 
X(NI,J,K)=XMAX 
ENDDO 

C --Bottom Surface 
K=2 
X(2,2,K)=XMIN 
Y(2,2,K)=YMIN 
DO I=3,NI-1 
X(I,2,K)=X(I-1,2,K)+DX(I-2) 
X(I,NJ,K)=X(I,2,K) 
Y(I,2,K)=YMIN 
Y(I,NJ,K)=YMAX 
ENDDO 

DO J=3,NJ-1 
Y(2,J,K)=Y(2,J-1,K)+DX(J-2) 
X(2,J,K)=XMIN 
Y(NI,J,K)=Y(2,J,K) 
X(NI,J,K)=XMAX 
ENDDO 

K=2 
DO I=2,NI 
DO J=2,NJ 
RI1=DBLE(I-2)/DBLE(NI-2) 
RI2=DBLE(NI-I)/DBLE(NI-2) 
X1(I,J,K)=RI1*X(NI,J,K)+RI2*X(2,J,K) 
Yl(I,J,K)=RI1*Y(NI,J,K)+RI2*Y(2,J,K) 
ENDDO 
ENDDO 

DO I=2,NI 
DO J=2,NJ 
RJ1=DBLE(J-2)/DBLE(NJ-2) 
RJ2=DBLE(NJ-J)/DBLE(NJ-2) 
X2(I,J,K)=RJ1*(X(I,NJ,K)-X1(I,NJ,K))+RJ2*(X(I,2,K)-Xl(I,2,K)) 
Y2(I,J,K)=RJ1*(Y(I,NJ,K)-Y1(I,NJ,K))+RJ2*(Y(I,2,K)-Yl(I,2,K)) 
ENDDO 
ENDDO 

DO I=2,NI 
DO J=2,NJ 
X(I,J,K)=X1(I,J,K)+X2(I,J,K) 
Y(I,J,K)=Y1(I,J,K)+Y2(I,J,K) 
ENDDO 
ENDDO 

Y(1,2,K)=Y(2,2,K)
X(1,2,K)=XMIN-DELTAX
X(NI+1,2,K)=Y(NI,2,K)
X(NI+1,2,K)=XMAX+DELTAX

K=2
DO J=2,NJ
DO I=2,NI
Z(I,J,K)=ZCURV(X(I,J,K),Y(I,J,K))
ENDDO 
ENDDO 

C --Top Surface 
K=NK
XO=XMIN
DO I=2,NI
XP=XO+DELTAX*(I-2)
X(I,2,K)=XP
Y(I,2,K)=YMIN
X(I,NJ,K)=XP
Y(I,NJ,K)=YMAX
ENDDO

DO I=2,NI 
YO = Y(I,2,K) 
DO J=3,NJ-1 
YP=YO+DELTAY 
YO=YP 
Y(I,J,K)=YP 
X(I,J,K)=X(I,NJ,K) 
X(1,J,K)=XMIN-DELTAX 
X(NIP1,J,K)=XMAX+DELTAX 
ENDDO 
ENDDO 

K=NK 
DO J=2, NJ 
DO I=2 ,NI 
Z(I,J,K)=ZMAX 
ENDDO 
ENDDO 

C --Front
I=NI
J=2
DO K=2,NK
Z(I,J,K)=ZMIN+DELTAZ*(K-2)
X(I,J,K)=XMAX 
Y(1,J,K)=YMIN 
ENDDO 

I=2 
J=2 
DO K=2,NK 
Z(I,J,K)=ZMIN+DELTAZ*(K-2) 
X(I,J,K)=XMIN 
Y(I,J,K)=YMIN 
ENDDO 

J=2
DO I=2,NI 
DO K=2,NK 
RI1=DBLE(I-2)/DBLE(NI-2) 
RI2=DBLE(NI-I)/DBLE(NI-2) 
X1(I,J,K)=RI1*(X(NI,J,K)+RI2*X(2,J,K)
Z1(I,J,K)=RI1*(Z(NI,J,K)+RI2*Z(2,J,K)
ENDDO 
ENDDO 

DO I=2,NI 
DO K=2,NK 
RJ1=DBLE(K-2)/DBLE(NK-2) 
RJ2=DBLE(NK-K)/DBLE(NK-2) 
X2(I,J,K)=RJ1*(X(I,J,NK)-X1(I,J,NK))+RJ2*(X(I,J,2)-X1(I,J,2)) 
Z2(I,J,K)=RJl*(Z(I,J,NK)-Z1(I,J,NK))+RJ2*(Z(I,J,2)-Z1(I,J,2)) 
ENDDO 
ENDDO 

DO I=2,NI
DO K=2,NK
X(I,J,K)=X1(I,J,K)+X2(I,J,K) 
Z(I,J,K)=Z1(I,J,K)+Z2(I,J,K) 
ENDDO 
ENDDO 

J=2 
DO K=2,NK 
DO I=2,NI 
        Y(I,J,K)=YMIN 
ENDDO 
ENDDO 

C --Back Surface 
I=NI 
J=NJ 
DO K=2,NK 
Z(I,J,K)=ZMIN+DELTAZ*(K-2) 
X(I,J,K)=XMAX 
Y(I,J,K)=YMAX 
ENDDO 

I=2 
J=NJ 
DO K=2,NK 
Z(I,J,K)=ZMIN+DELTAZ*(K-2) 
X(I,J,K)=XMIN 
Y(I,J,K)=YMAX 
ENDDO 

J=NJ 
DO I=2,NI 
DO K=2,NK 
RI1=DBLE(I-2)/DBLE(NI-2) 
RI2=DBLE(NI-I)/DBLE(NI-2) 
X1(I,J,K)=RI1*X(NI,J,K)+RI2*X(2,J,K) 
Z1(I,J,K)=RI1*Z(NI,J,K)+RI2*Z(2,J,K) 
ENDDO 
ENDDO 

DO I=2 ,NI 
DO K=2,NK 
RJ1=DBLE(K-2)/DBLE(NK-2) 
RJ2=DBLE(NK-K)/DBLE(NK-2) 
X2(I,J,K)=RJ1*(X(I,J,NK)-X1(I,J,NK))+RJ2*(X(I,J,2)-X1(I,J,2))
Z2(I,J,K)=RJl*(Z(I,J,NK)-Zl(I,J,NK))+RJ2*(Z(I,J,2)-Z1(I,J,2)) 
ENDDO 
ENDDO 

DO I=2,NI 
DO K=2,NK 
X(I,J,K)=X1(I,J,K)+X2(I,J,K) 
Z(I,J,K)=Z1(I,J,K)+Z2(I,J,K) 
ENDDO 
ENDDO 

J=NJ
DO K=2,NK
DO I=2,NI
        Y(I,J,K)=YMAX 
ENDDO 
ENDDO 

C --Left Surface 
I=2 
DO J=2 ,NJ 
DO K=2,NK 
RI1=DBLE(J-2)/DBLE(NJ-2) 
RI2=DBLE(NJ-J)/DBLE(NJ-2) 
Y1(I,J,K)=RI1*Y(I,NJ,K)+RI2*Y(I,2,K) 
Z1(I,J,K)=RI1*Z(I,NJ,K)+RI2*Z(I,2,K) 
ENDDO 
ENDDO 

DO J=2,NJ 
DO K=2,NK 
RJ1=DBLE(K-2)/DBLE(NK-2) 
RJ2=DBLE(NK-K)/DBLE(NK-2) 
Y2(I,J,K)=RJ1*(Y(I,J,NK)-Y1(I,J,NK))+RJ2*(Y(I,J,2)-Y1(I,J,2))
Z2(I,J,K)=RJ1*(Z(I,J,NK)-Z1(I,J,NK))+RJ2*(Z(I,J,2)-Z1(I,J,2)) 
ENDDO
ENDDO

DO J=2,NJ 
DO K=2,KK 
Y(I,J,K)=Y1(I,J,K)+Y2(I,J,K) 
Z(I,J,K)=Z1(I,J,K)+Z2(I,J,K) 
ENDDO 
ENDDO 

I=2 
DO K=2, NK 
DO J=2, NJ 
X(I,J,K)=XMIN 
ENDDO 
ENDDO 

C --Right Surface 
I=NI
DO J=2,NJ
DO K=2,NK
RI1=DBLE(J-2)/DBLE(NJ-2) 
RI2=DBLE(NJ-J)/DBLE(NJ-2) 
Y1(I,J,K)=RI1*Y(I,NJ,K)+RI2*Y(I,2,K) 
Z1(I,J,K)=RI1*Z(I,NJ,K)+RI2*Z(I,2,K)
ENDDO 
ENDDO 

DO J=2 ,NJ 
DO K=2,NK 
RJ1=DBLE(K-2)/DBLE(NK-2) 
RJ2=DBLE(NK-K)/DBLE(NK-2) 
Y2(I,J,K)=RJ1*(Y(I,J,NK)-Y1(I,J,NK))+RJ2*(Y(I,J,2)-Y1(I,J,2)) 
Z2(I,J,K)=RJ1*(Z(I,J,NK)-Z1(I,J,NK))+RJ2*(Z(I,J,2)-Z1(I,J,2)) 
ENDDO 
ENDDO 

DO J=2,NJ 
DO K=2,NK 
Y(I,J,K)=Y1(I,J,K)+Y2(I,J,K) 
Z(I,J,K)=Z1(I,J,K)+Z2(I,J,K) 
ENDDO 
ENDDO 

I=NI 
DO K=2,NK 
DO J=2,NJ 
X(I,J,K)=XMAX 
ENDDO 
ENDDO 

C --Initial Grid

DO I=3,NI-1
DO J=3,NJ-1 
        DELTAZ=(Z(I,J,NK)-Z(I,J,2))/DBLE(NK-2)
DO K=3,NK-1
Z(I,J,K)=Z(I,J,K-1)+DELTAZ
X(I,J,K)=X(I,J,2) 
Y(I,J,K)=Y(I,J,2) 
ENDDO 
ENDDO 
ENDDO 

C --Save the grid to PlotMtv format file 
CALL GRDEXT(X,Y,Z)
CALL SAVEGRD(X,Y,Z)
CALL SAVEMTV(X,Y,Z) 

RETURN 
END 


SUBROUTINE GRDEXT (X,Y,Z) 
include "common.h" 
REAL*8 X(NIP1,NIP1,NKP1) 
REAL*8 Y(NIP1,NIP1,NKP1) 
REAL*8 Z(NIP1,NIP1,NKP1) 
REAL*8 tmp1 
INTEGER I, J ,K 

DO K=2,NK
DO J=2,NJ
X(1,J,K)=2.0D0*X(2,J,K)-X(3,J,K) 
X(NIP1,J,K)=2.0D0*X(NI,J,K)-X(NI-1,J,K)
Y(1,J,K)=Y(3,J,K)
Y(NIP1,J,K)=Y(NI-1,J,K) 
Z(I,J,K)=Z(3,J,K) 
Z(NIP1,J,K)=Z(NI-1,J,K) 
ENDDO 
ENDDO 

DO K=2, NK 
DO I=2 ,NI 
X(I,1,K)=X(I,3,K) 
X(I,NIP1,K)=X(I,NJ-1,K) 
Y(I,1,K)=2.0D0*Y(I,2,K)-Y(I,3,K) 
Y(I,NIP1,K)=2.0D0*Y(I,NJ,K)-Y(I,NJ-1,K) 
Z(I,1,K)=Z(I,3,K) 
Z(I,NJP1,K)=Z(I,NJ-1,K) 
ENDDO 
ENDDO 

C ---Corner Points 

DO K=2 ,NK 
tmp1=2.0D0*X(2,1,K)-X(3,1,K) 
X(1,1,K)=2.0D0*X(2,2,K)-X(3,3,K) 
tmp1=2.0D0*X(2,NJP1,K)-X(3,NJP1,K) 
X(1,NJP1,K)=2.0D0*X(2,NJ,K)-X(3,NJ-1,K) 
tmp1=2.0D0*X(NI,1,K)-X(NI-1,1,K) 
X(NIP1,1,K)=2.0D0*X(NI,2,K)-X(NI-1,3,K) 
tmp1 = 2.0D0*X(NI,NJP1,K)-X(NI-1,NJP1,K) 
X(NIP1,NIP1,K)=2.0DO*X(NI,NJ,K)-X(NI-l,NJ-1,K) 

tmp1 = 2.0D0*Y(1,2,K)-Y(1,3,K) 
Y(1,1,K)= 0.5D0*(tmp1+0.5D0*(Y(2,1,K)+Y(3,1,K))) 
tmp1 = 2.0D0*Y(1,NI,K)-Y(1,NI-1,K) 
Y(1,NJP1,K)= 0.5D0*(tmp1+0.5D0*(Y(2,NJP1,K)+Y(3,NJP1,K))) 
tmp1 = 2.0D0*Y(NIP1,2,K)-Y(NIP1,3,K) 
Y(NIP1,1,K)= 0.5D0*(tmp1+0.5D0*(Y(NI,1,K)+Y(NI-1,1,K))) 
tmp1 = 2.0D0*Y(NIP1,NI,K)-Y (NIP1,NI-1,K) 
Y(NIP1,NJP1,K)=0.5D0*(tmp1+0.5D0*(Y(NI,NJP1,K)+Y(NI-1,NJP1,K))) 

Z(1,1,K)=Z(3,3,K) 
Z(1,NJP1,K)=Z(3,NJ-1,K) 
Z(NIP1,1,K)=Z(NI-1,3,K) 
Z(NIP1,NJP1,K)=Z(NI-1,NJ-1,K) 
ENDDO 
RETURN 
END 

SUBROUTINE SEESPC(X,Y,Z)
include "common.h"
REAL*8 X(NIP1,NJP1,NKP1),Y(NIP1,NJP1,NKP1),Z(NIP1,NJP1,NKP1)
REAL*8 ADS,DX,DY,DZ,MAXADS
INTEGER I,J
MAXADS=0.0D0
DO J=2,NJ
DO I=2,NI
DX=X(I,J,3)-X(I,J,2)
DY=Y(I,J,3)-Y(I,J,2) 
DZ=Z(I,J,3)-Z(I,J,2) 
ADS=DSQRT(DX*DX+DY*DY+DZ*DZ) 
MAXADS=MAX(MAXADS,ADS) 
ENDDO 
ENDDO 
PRINT*,'MAX DS  ',MAXADS 
RETURN 
END

SUBROUTINE POW(X,N,PW) 
REAL*8 X,TMP,PW 
INTEGER N 
TMP=N*LOG(X) 
PW=EXP(TMP) 
RETURN 
END 

SUBROUTINE SAVEMTV( X,Y, Z) 
include "common.h" 
REAL*8 X(NIP1,NIP1,NKP1) 
REAL*8 Y (NIP1,NIP1,NKP1) 
REAL*8 Z (NIP1,NIP1,NKP1) 
INTEGER I, J, K 

OPEN(23,FILE='3dgrid.mtv') 
WRITE(23,998)'$ DATA = CURVE3D' 
WRITE(23,998)'% axisscale = False' 
WRITE(23,998)'% axislabel = False' 

I=16
DO K=2,NK
DO J=2,NJ
WRITE(23,*),X(I,J,K),' ',Y(I,J,K),' ',Z(I,J,K)
ENDDO 
WRITE(23,*) 
ENDDO 
WRITE(23,*) 

DO J=2 ,NJ 
DO K=2,NK 
WRITE(23,*),X(I,J,K),' ',Y(I,J,K),' ',Z(I,J,K) 
ENDDO 
WRITE(23,*) 
ENDDO 
WRITE(23,*) 

J=16 
DO K=2,NK 
DO I=2,NI 
WRITE(23,*),X(I,J,K),' ',Y(I,J,K),' ',Z(I,J,K)
ENDDO 
WRITE(23,*) 
ENDDO 
WRITE(23.*) 

K=2
DO J=12,NJ
DO I=16,NI
WRITE(23,*),X(I,J,K),' ',Y(I,J,K),' ',Z(I,J,K)
ENDDO
WRITE(23,*) 
ENDDO 
WRITE(23,*) 

DO I=16,NI
DO J=2,NJ
WRITE(23,*),X(I,J,K),' ',Y(I,J,K),' ',Z(I,J,K) 
ENDDO 
WRITE(23,*) 
ENDDO 
WRITE(23,*) 

CLOSE (23) 
998 	FORMAT(A) 
RETURN 
END 

SUBROUTINE SAVEGRD(X,Y,Z) 
include "common.h" 
REAL*8 X(NIP1,NJP1,NKP1),Y(NIP1,NIP1,NKP1),Z(NIP1,NIP1,NKP1) 
INTEGER I,J,K 
OPEN(24,FILE='3dgrid.plt') 
WRITE(24,*)'VARIABLES = "X","Y","Z","T"'
WRITE(24,*)'ZONE F=POINT,I=',NI-1,',J=',NJ-1,',K=',NK-1 
DO K=2,NK 
DO J=2,NJ 
DO I=2,NI 
        WRITE(24,*),X(I,J,K),' ',Y(I,J,K),' ',Z(I,J,K),'0.0'
ENDDO 
ENDDO 
ENDDO 
WRITE(24,*) 
CLOSE(24) 
1001 	FORMAT(A) 
RETURN 
END 




